GASPARONI dataset

1 Data retrival

library(dplyr)
library(minfi)

## Warning: package 'S4Vectors' was built under R version 3.6.3

## Warning: package 'GenomeInfoDb' was built under R version 3.6.3

## Warning: package 'DelayedArray' was built under R version 3.6.3

library(IlluminaHumanMethylation450kmanifest)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)

cohort <- "GASPARONI"
data.dir <- file.path("../DATASETS",cohort) 
data.dir.table <- "../DATASETS/Summary_Table/" 
data.dir.raw <- "../../coMethDMR_metaAnalysis/code_validation/Meta_analysis_code/DATASETS/GASPARONI/step2_read_minfi/"
data.dir.bsfilter <- file.path(data.dir,"step2_bsConvFilter/") 
data.dir.clinical.filter <- file.path(data.dir,"step3_clinical_available_filtering/") 
data.dir.probes.qc <- file.path(data.dir,"step4_probesQC_filtering/") 
data.dir.probes.normalization <- file.path(data.dir,"step5_normalization/") 
data.dir.pca <- file.path(data.dir,"step6_pca_filtering/") 
data.dir.neuron <- file.path(data.dir,"step7_neuron_comp/") 
data.dir.single.cpg.pval <- file.path(data.dir,"step8_single_cpg_pval/") 
data.dir.residuals <- file.path(data.dir,"step9_residuals/") 
data.dir.median <- file.path(data.dir,"step10_median/") 
for(p in grep("dir",ls(),value = T)) dir.create(get(p),recursive = TRUE,showWarnings = FALSE)

2 Confirm sex status sing DNA methylation values

# Create a MethylSet object from RGSet
MSet <- preprocessRaw(RGSet)

# Create [Genomic]MethylSet object
GMset <- mapToGenome(MSet)

# Get predicted sex status
estSex <- getSex(GMset)

# Compare predicted gender with phenotype gender
realSex <- data.frame(
  sample = paste(phenoData$geo_accession,
                 phenoData$sentrix_id.ch1,
                 phenoData$sentrix_position.ch1,
                 sep = "_"),
  realSex = phenoData$Sex.ch1,
  stringsAsFactors = FALSE
)
compareSex <- merge(
  as.data.frame(estSex), realSex,
  by.x = "row.names",
  by.y = "sample")

identical(compareSex$predictedSex, compareSex$realSex)

## [1] TRUE

print(compareSex)

##                       Row.names     xMed     yMed predictedSex realSex
## 1  GSM2808939_5854945043_R03C01 13.39486 10.71210            F       F
## 2  GSM2808940_5854945043_R04C01 12.63866 12.88729            M       M
## 3  GSM2808941_5854945043_R05C01 13.04055 10.40939            F       F
## 4  GSM2808942_5854945043_R06C01 13.28677 10.65374            F       F
## 5  GSM2808943_5854945043_R01C02 12.39178 12.63091            M       M
## 6  GSM2808944_5854945043_R02C02 12.60409 12.81878            M       M
## 7  GSM2808945_5854945043_R03C02 13.41924 10.62388            F       F
## 8  GSM2808946_5854945043_R04C02 13.39158 10.54978            F       F
## 9  GSM2808947_5854945043_R05C02 13.24926 10.63209            F       F
## 10 GSM2808948_5854945043_R06C02 12.39044 12.67881            M       M
## 11 GSM2808949_5854945045_R01C01 12.42285 12.74051            M       M
## 12 GSM2808950_5854945045_R02C01 12.37368 12.66977            M       M
## 13 GSM2808951_5854945045_R03C01 12.73450 13.00885            M       M
## 14 GSM2808952_5854945045_R04C01 13.26298 10.47978            F       F
## 15 GSM2808953_5854945045_R05C01 13.26408 10.53236            F       F
## 16 GSM2808954_5854945045_R06C01 13.08688 10.48382            F       F
## 17 GSM2808955_5854945045_R01C02 12.11309 12.46097            M       M
## 18 GSM2808956_5854945045_R02C02 12.11667 12.45635            M       M
## 19 GSM2808957_5854945045_R03C02 12.65486 12.87517            M       M
## 20 GSM2808958_5854945045_R04C02 12.56332 12.83624            M       M
## 21 GSM2808959_5854945045_R05C02 12.48066 12.75666            M       M
## 22 GSM2808961_5854945056_R01C01 13.13627 10.39339            F       F
## 23 GSM2808962_5854945056_R02C01 12.49035 12.75822            M       M
## 24 GSM2808963_5854945056_R03C01 13.35996 10.64926            F       F
## 25 GSM2808964_5854945056_R04C01 12.71928 13.02574            M       M
## 26 GSM2808965_5854945056_R05C01 12.85214 13.06623            M       M
## 27 GSM2808966_5854945056_R06C01 13.42528 10.78259            F       F
## 28 GSM2808967_5854945056_R01C02 12.45108 12.68650            M       M
## 29 GSM2808968_5854945056_R02C02 12.61102 12.85428            M       M
## 30 GSM2808969_5854945056_R03C02 12.77345 13.07364            M       M
## 31 GSM2808970_5854945056_R04C02 13.47763 10.53625            F       F
## 32 GSM2808971_5854945056_R05C02 13.15014 10.55651            F       F
## 33 GSM2808972_5854945056_R06C02 12.49935 12.75488            M       M
## 34 GSM2808973_5854945057_R01C01 13.30592 10.39714            F       F
## 35 GSM2808974_5854945057_R02C01 12.50767 10.04029            F       F
## 36 GSM2808975_5854945057_R03C01 13.18627 10.48432            F       F
## 37 GSM2808976_5854945057_R04C01 13.44016 10.60918            F       F
## 38 GSM2808977_5854945057_R05C01 12.77582 13.02998            M       M
## 39 GSM2808978_5854945057_R06C01 13.38262 10.77602            F       F
## 40 GSM2808979_5854945057_R01C02 13.23937 10.24317            F       F
## 41 GSM2808980_5854945057_R02C02 13.23691 10.37286            F       F
## 42 GSM2808981_5854945057_R03C02 12.91663 10.19660            F       F
## 43 GSM2808982_5854945057_R05C02 13.40886 10.59898            F       F
## 44 GSM2808983_5854945057_R06C02 13.40221 10.63617            F       F
## 45 GSM2808984_5854945010_R01C01 12.04644 12.24510            M       M
## 46 GSM2808985_5854945010_R02C01 12.94782 10.26562            F       F
## 47 GSM2808986_5854945010_R03C01 13.25820 10.25562            F       F
## 48 GSM2808987_5854945010_R06C01 12.80886 13.01280            M       M
## 49 GSM2808988_5854945010_R01C02 12.31217 12.51607            M       M
## 50 GSM2808989_5854945010_R02C02 13.18936 10.67684            F       F
## 51 GSM2808990_5854945010_R03C02 12.75937 12.95074            M       M
## 52 GSM2808991_5854945010_R04C02 12.88741 13.10648            M       M
## 53 GSM2808992_5854945010_R05C02 12.52613 12.74398            M       M
## 54 GSM2808993_5854945010_R06C02 13.42423 11.03411            F       F
## 55 GSM2808994_5854945021_R01C01 12.96208 10.15924            F       F
## 56 GSM2808995_5854945021_R02C01 12.33330 12.53453            M       M
## 57 GSM2808996_5854945021_R03C01 13.20823 10.48684            F       F
## 58 GSM2808997_5854945021_R04C01 12.59397 12.82137            M       M
## 59 GSM2808998_5854945021_R06C01 13.45950 10.85526            F       F
## 60 GSM2808999_5854945021_R01C02 12.80413 10.25680            F       F

3 Data QC

3.1 Bilsufite conversion filtering

Removing samples with bisulfiteConversion lower than 88.

library(wateRmelon)
library(RPMM)
#library(sesame)
#library(sesameData)

phenoData <- phenoData[match(substr(colnames(RGSet),1,10), phenoData$geo_accession),]
nb.samples <- nrow(phenoData)
nb.female.samples <- sum(phenoData$Sex.ch1 == "F")
nb.male.samples <- sum(phenoData$Sex.ch1 == "M")

bs <- data.frame(bisulfiteConversion = bscon(RGSet))
bsFilteredOut <- row.names(bs)[bs$bisulfiteConversion < 88]
RGSet <- RGSet[,!colnames(RGSet) %in% bsFilteredOut]
phenoData <- phenoData[match(substr(colnames(RGSet),1,10), phenoData$geo_accession),]

nb.samples.bc.filtered <- nrow(phenoData)
nb.female.samples.bc.filtered <- sum(phenoData$Sex.ch1 == "F")
nb.male.samples.bc.filtered <- sum(phenoData$Sex.ch1 == "M")

ggpubr::gghistogram(bs$bisulfiteConversion,xlab = "bisulfite Conversion" )

## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.

save(RGSet,
     nb.samples,
     nb.female.samples,
     nb.male.samples,
     bs,
     phenoData,
     nb.samples.bc.filtered,
     nb.female.samples.bc.filtered,
     nb.male.samples.bc.filtered,
     file = paste0(data.dir.bsfilter, "/RGSet_bsfiltered.rda"))

3.2 Clinical data filtering

phenoData$braak_stage.ch1 <- phenoData$braak_stage.ch1 %>% as.numeric()
phenoData$age.ch1 <- phenoData$age.ch1 %>% as.numeric()
### Subset rows and columns
pheno_df <- phenoData  %>% as.data.frame() %>%
  dplyr::filter(
    source_name_ch1 == "Frontal Cortex" &
      !is.na(phenoData$braak_stage.ch1) &
      phenoData$characteristics_ch1 == "cell type: bulk"
  ) %>% dplyr::select(
    c(
      "geo_accession",
      "donor_id.ch1",
      "sentrix_id.ch1",
      "age.ch1",
      "Sex.ch1",
      "braak_stage.ch1"
    )
  )

### Rename vars
colnames(pheno_df) <- c(
  "sample", "subject.id", "slide", "age.brain", "sex", "stage"
)
nb.samples.with.clinical <- nrow(pheno_df)
nb.female.samples.with.clinical <- sum(pheno_df$sex == "F")
nb.male.samples.with.clinical <- sum(pheno_df$sex == "M")

## phenotype dataset
save(RGSet,
     nb.samples.with.clinical,
     nb.female.samples.with.clinical,
     nb.male.samples.with.clinical,
     pheno_df,
     file = paste0(data.dir.clinical.filter, "/gasparoni_bs_and_clinical_filtered.rda"))

3.3 Probes QC

Input:

• RGSet.RDS
• beta_mat.RDS

Output:

• beta_CG_XY_SNPfiltered.RDS

library(minfi)
library(DMRcate)

## Warning: replacing previous import 'minfi::getMeth' by 'bsseq::getMeth' when
## loading 'DMRcate'

library(sesame)

### Find which chromosome each probe is on
beta_mat <- beta_mat <- getBeta(RGSet)

probes.info <- sesameDataGet("HM450.hg19.manifest")
probes.info <- probes.info[row.names(beta_mat) %>% as.character()] %>%
  as.data.frame %>%
  dplyr::select(c("seqnames","start","end"))
probes.info$seqnames <- as.character(probes.info$seqnames)

nb.probes <- nrow(probes.info)
nb.chrAuto.probes <- sum(probes.info$seqnames %in% paste0("chr", 1:22))
nb.chrX.probes <- sum(probes.info$seqnames == "chrX")
nb.chrY.probes <- sum(probes.info$seqnames == "chrY")
nb.chrM.probes <- sum(probes.info$seqnames == "chrM")

### subset to probes with detection P <= 0.01
detP <- detectionP(RGSet, type = "m+u")
failed.01 <- detP > 0.01
passedProbes <- rownames(failed.01)[rowMeans(failed.01) == 0] 
beta_mat <- beta_mat[passedProbes, ]
probes.info <- probes.info[row.names(probes.info) %in% row.names(beta_mat),]

nb.probes.detectP <- nrow(probes.info)
nb.chrAuto.probes.detectP <- sum(probes.info$seqnames %in% paste0("chr", 1:22))
nb.chrX.probes.detectP <- sum(probes.info$seqnames == "chrX")
nb.chrY.probes.detectP <- sum(probes.info$seqnames == "chrY")
nb.chrM.probes.detectP <- sum(probes.info$seqnames == "chrM")

### keep only probes that start with "cg"
beta_mat <- beta_mat[grep("cg",rownames(beta_mat)),]
probes.info <- probes.info[row.names(probes.info) %in% row.names(beta_mat),]

nb.probes.detectP.cg <- nrow(probes.info)
nb.chrAuto.probes.detectP.cg <- sum(probes.info$seqnames %in% paste0("chr", 1:22))
nb.chrX.probes.detectP.cg <- sum(probes.info$seqnames == "chrX")
nb.chrY.probes.detectP.cg <- sum(probes.info$seqnames == "chrY")
nb.chrM.probes.detectP.cg <- sum(probes.info$seqnames == "chrM")

### drop probes where SNP with MAF >= 0.01 in the last 5 bp of the probe
beta_mat <- rmSNPandCH(
  object = beta_mat,
  dist = 5, 
  mafcut = 0.01, 
  and = TRUE,
  rmcrosshyb = FALSE,
  rmXY = FALSE
)
probes.info <- probes.info[row.names(probes.info) %in% row.names(beta_mat),]

nb.probes.cg.dmrcate <- nrow(probes.info)
nb.chrAuto.probes.cg.dmrcate <- sum(probes.info$seqnames %in% paste0("chr", 1:22))
nb.chrX.probes.cg.dmrcate <- sum(probes.info$seqnames == "chrX")
nb.chrY.probes.cg.dmrcate <- sum(probes.info$seqnames == "chrY")
nb.chrM.probes.cg.dmrcate <- sum(probes.info$seqnames == "chrM")

### drop probes in chrM
probes.info <- probes.info[probes.info$seqnames != "chrM",]
beta_mat <- beta_mat[
  row.names(beta_mat) %in% row.names(probes.info),
]

nb.probes.dmrcate.chrM <- nrow(probes.info)
nb.chrAuto.probes.dmrcate.chrM <- sum(probes.info$seqnames %in% paste0("chr", 1:22))
nb.chrX.probes.dmrcate.chrM <- sum(probes.info$seqnames == "chrX")
nb.chrY.probes.dmrcate.chrM <- sum(probes.info$seqnames == "chrY")

save(
  pheno_df,
  beta_mat,
  nb.probes.detectP,
  nb.chrAuto.probes.detectP,
  nb.chrX.probes.detectP,
  nb.chrY.probes.detectP,
  nb.chrM.probes.detectP,
  nb.probes.detectP.cg,
  nb.chrAuto.probes.detectP.cg,
  nb.chrX.probes.detectP.cg,
  nb.chrY.probes.detectP.cg,
  nb.chrM.probes.detectP.cg,
  nb.probes.cg.dmrcate,
  nb.chrAuto.probes.cg.dmrcate,
  nb.chrX.probes.cg.dmrcate,
  nb.chrY.probes.cg.dmrcate,
  nb.chrM.probes.cg.dmrcate,
  nb.probes.dmrcate.chrM,
  nb.chrAuto.probes.dmrcate.chrM,
  nb.chrX.probes.dmrcate.chrM,
  nb.chrY.probes.dmrcate.chrM,
  file = paste0(data.dir.probes.qc, "/beta_CG_XY_SNPfiltered_mat.rda")
)

4 Normalization

• Quantile normalization and BMIQ normalization

Input:

• beta_CG_XY_SNPfiltered_mat.RDS
• RGSet.RDS
• pheno_df.RDS
• full.annot.RDS

Output:

• bs.csv
• pheno_df.RDS
• QNBMIQ.RDS

4.1 Quantile normalization

4.1.1 Check if sex chrom have higher/lower intensities than autosomes

4.1.1.1 Create function to extract separate matrix for each chromosome

library(sesame)
library(ggplot2)

probes.info <- sesameDataGet("HM450.hg19.manifest")
probes.info <- probes.info[row.names(beta_mat) %>% as.character()] %>%
  as.data.frame %>%
  dplyr::select(c("seqnames","start","end"))
probes.info$seqnames <- as.character(probes.info$seqnames)

findBetaChr <- function(data, chrom){
  data %>%
    as.data.frame() %>% # has to turn matrix into df for next step
    tibble::rownames_to_column() %>% # turn rownames to a column, so rownames won't be deleted after filtering rows in next step
    filter(rowname %in% row.names(probes.info[probes.info$seqnames == chrom,])) %>%
    tibble::column_to_rownames() %>%
    as.matrix()
}

4.1.1.2 Compare beta values between gender

4.1.1.2.1 on chrX

betaChrX <- findBetaChr(data = beta_mat, chrom = "chrX")
betaChrX_long <- data.frame(
  beta = as.vector(betaChrX),
  sample = rep(substr(colnames(betaChrX), 1,10), each = nrow(betaChrX)),
  stringsAsFactors = FALSE
)
betaChrX_long <- merge(
  betaChrX_long, pheno_df[, c("sample", "sex")],
  by = "sample",
  sort = FALSE
)
  
ggplot(betaChrX_long,
       aes(x = sample, y = beta, fill = sex)) +
  stat_boxplot(geom ='errorbar', width = 1, linetype = 1) +
  geom_boxplot(width = 1, alpha = 1, outlier.shape = 1, outlier.size = 2) +
  scale_fill_grey(start=1, end=0.6) +
  labs(x = "sex", y = "DNA methylation beta values",
       title = "DNA methylation level by sex on chromosome X") +
  theme_bw() +
  theme(legend.position = "none", axis.text.x = element_blank()) +
  facet_wrap(~sex)

4.1.1.2.2 on chrY

betaChrY <- findBetaChr(data = beta_mat, chrom = "chrY")
betaChrY_long <- data.frame(
  beta = as.vector(betaChrY),
  sample = rep(substr(colnames(betaChrY), 1,10), each = nrow(betaChrY)),
  stringsAsFactors = FALSE
)
betaChrY_long <- merge(
  betaChrY_long, pheno_df[, c("sample", "sex")],
  by = "sample",
  sort = FALSE
)

ggplot(betaChrY_long,
       aes(x = sample, y = beta, fill = sex)) +
  stat_boxplot(geom ='errorbar', width = 1, linetype = 1) +
  geom_boxplot(width = 1, alpha = 1, outlier.shape = 1, outlier.size = 2) +
  scale_fill_grey(start=1, end=0.6) +
  labs(x = "sex", y = "DNA methylation beta values",
       title = "DNA methylation level by sex on chromosome Y") +
  theme_bw() +
  theme(legend.position = "none", axis.text.x = element_blank()) +
  facet_wrap(~sex)

4.1.1.3 Boxplot of methylation by chromosomes

4.1.1.3.1 overall

chrAuto <- paste0("chr", 1:22)
betaChrAuto_ls <- lapply(seq_along(chrAuto), function(i){
  dat <- findBetaChr(data = beta_mat, chrom = chrAuto[i])
  data.frame(beta = as.vector(dat), chrom = "autosomes", stringsAsFactors = FALSE)
})
betaChrAuto_df <- do.call(rbind, betaChrAuto_ls)

chrSex <- paste0("chr", c("X", "Y"))
betaChrSex_ls <- lapply(seq_along(chrSex), function(i){
  dat <- findBetaChr(data = beta_mat, chrom = chrSex[i])
  data.frame(beta = as.vector(dat), chrom = chrSex[i], stringsAsFactors = FALSE)
})
betaChrSex_df <- do.call(rbind, betaChrSex_ls)

betaChr_df <- rbind(betaChrAuto_df, betaChrSex_df)

# # No enough memory to load all the data points and plot the figure, so use function boxplot instead
# ggplot(betaChr_df,
#        aes(x = chrom, y = beta, fill = chrom)) +
#   stat_boxplot(geom ='errorbar', width = 0.4, linetype = 1) +
#   geom_boxplot(width = 0.4, alpha = 1, outlier.shape = 1, outlier.size = 2) +
#   scale_fill_grey(start=1, end=0.6) +
#   labs(x = "chromosome type", y = "DNA methylation beta values",
#        title = "DNA methylation level vs. chromosome type") +
#   theme_bw()

boxplot(
  beta ~ chrom, data = betaChr_df, ylab = "DNA methylation beta values",
  main = "DNA methylation on chromosomes")

pdf(file = paste0(data.dir.probes.normalization, "/beta_vs_chrom.pdf"))
boxplot(
  beta ~ chrom, data = betaChr_df, ylab = "DNA methylation beta values",
  main = "DNA methylation on chromosomes")
dev.off()

## quartz_off_screen 
##                 2

4.1.1.3.2 by gender

### split matrix by sex first
beta_mat_female <- beta_mat[ 
  ,substr(colnames(beta_mat), 1, 10) %in% pheno_df$sample[pheno_df$sex == "F"]]
beta_mat_male <- beta_mat[ 
  ,substr(colnames(beta_mat), 1, 10) %in% pheno_df$sample[pheno_df$sex == "M"]]

### for females
chrAuto <- paste0("chr", 1:22)
female_auto_ls <- lapply(seq_along(chrAuto), function(i){
  dat <- findBetaChr(data = beta_mat_female, chrom = chrAuto[i])
  data.frame(beta = as.vector(dat), chrom = "autosomes", stringsAsFactors = FALSE)
})
female_auto_df <- do.call(rbind, female_auto_ls)

female_sex_df <- findBetaChr(data = beta_mat_female, chrom = "chrX")
female_sex_df <- data.frame(beta = as.vector(female_sex_df), chrom = "chromosome X", stringsAsFactors = FALSE)

female_df <- rbind(female_auto_df, female_sex_df)

boxplot(
  beta ~ chrom, data = female_df, ylab = "DNA methylation beta values",
  main = "DNA methylation on chromosomes for females")

### for males
chrAuto <- paste0("chr", 1:22)
male_auto_ls <- lapply(seq_along(chrAuto), function(i){
  dat <- findBetaChr(data = beta_mat_male, chrom = chrAuto[i])
  data.frame(beta = as.vector(dat), chrom = "autosomes", stringsAsFactors = FALSE)
})
male_auto_df <- do.call(rbind, male_auto_ls)

male_X_df <- findBetaChr(data = beta_mat_male, chrom = "chrX")
male_X_df <- data.frame(beta = as.vector(male_X_df), chrom = "chromosome X", stringsAsFactors = FALSE)

male_Y_df <- findBetaChr(data = beta_mat_male, chrom = "chrY")
male_Y_df <- data.frame(beta = as.vector(male_Y_df), chrom = "chromosome Y", stringsAsFactors = FALSE)

male_df <- rbind(male_auto_df, male_X_df, male_Y_df)

boxplot(
  beta ~ chrom, data = male_df, ylab = "DNA methylation beta values",
  main = "DNA methylation on chromosomes for males")

### save plots
pdf(file = paste0(data.dir.probes.normalization, "/beta_vs_chrom_by_gender.pdf"))

boxplot(
  beta ~ chrom, data = female_df, ylab = "DNA methylation beta values",
  main = "DNA methylation on chromosomes for females")

boxplot(
  beta ~ chrom, data = male_df, ylab = "DNA methylation beta values",
  main = "DNA methylation on chromosomes for males")

dev.off()

## quartz_off_screen 
##                 2

4.1.1.4 Boxplot of methylation by gender on autosomes, X, Y

library(quantro)

4.1.1.4.1 on autosomes

chrAuto <- paste0("chr", 1:22)
betaChrAuto_ls <- lapply(seq_along(chrAuto), function(i){findBetaChr(data = beta_mat, chrom = chrAuto[i])})
betaChrAuto_df <- do.call(rbind, betaChrAuto_ls)
matboxplot(betaChrAuto_df,
           groupFactor = pheno_df$sex,
           xaxt = "n",
           main = "Beta Values (autosomes) - before normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

4.1.1.4.2 on chrX

betaChrX_df <- findBetaChr(data = beta_mat, chrom = "chrX")
matboxplot(betaChrX_df,
           groupFactor = pheno_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome X) - before normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

4.1.1.4.3 on chrY

betaChrY_df <- findBetaChr(data = beta_mat,chrom = "chrY")
matboxplot(betaChrY_df,
           groupFactor = pheno_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome Y) - before normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

4.1.1.4.4 save plots

pdf (
  paste0(data.dir.probes.normalization, "/boxPlotBeforeNormalization.pdf")
)

matboxplot(betaChrAuto_df,
           groupFactor = pheno_df$sex,
           xaxt = "n",
           main = "Beta Values (autosomes) - before normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

matboxplot(betaChrX_df,
           groupFactor = pheno_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome X) - before normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

matboxplot(betaChrY_df,
           groupFactor = pheno_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome Y) - before normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

dev.off()

## quartz_off_screen 
##                 2

4.1.2 normalization by group

library(lumi)

4.1.2.1 split beta_mat into sub-matrics by females autosomes, females X, males autosomes, males X, males Y

### split matrix by sex first
beta_mat_female <- beta_mat[ 
  ,substr(colnames(beta_mat), 1, 10) %in% pheno_df$sample[pheno_df$sex == "F"]]
beta_mat_male <- beta_mat[ 
  ,substr(colnames(beta_mat), 1, 10) %in% pheno_df$sample[pheno_df$sex == "M"]]

### split female beta matrix by chromosome (auto, X)
chrAuto <- paste0("chr", 1:22)
beta_mat_female_auto_ls <- lapply(seq_along(chrAuto), function(i){
  findBetaChr(data = beta_mat_female, chrom = chrAuto[i])})
beta_mat_female_auto <- do.call(rbind, beta_mat_female_auto_ls)

beta_mat_female_X <- findBetaChr(data = beta_mat_female, chrom = "chrX")

### split male beta matrix by chromosome (auto, X, Y)
chrAuto <- paste0("chr", 1:22)
beta_mat_male_auto_ls <- lapply(seq_along(chrAuto), function(i){
  findBetaChr(data = beta_mat_male, chrom = chrAuto[i])})
beta_mat_male_auto <- do.call(rbind, beta_mat_male_auto_ls)

beta_mat_male_X <- findBetaChr(data = beta_mat_male, chrom = "chrX")

beta_mat_male_Y <- findBetaChr(data = beta_mat_male, chrom = "chrY")

4.1.2.2 normalization

### female autosomes
betaQN_female_auto <- lumiN(x.lumi = beta_mat_female_auto, method = "quantile")

## Perform quantile normalization ...

dim(betaQN_female_auto)

## [1] 446784     30

### female chromosome X
betaQN_female_X <- lumiN(x.lumi = beta_mat_female_X, method = "quantile")

## Perform quantile normalization ...

dim(betaQN_female_X)

## [1] 10618    30

### male autosomes
betaQN_male_auto <- lumiN(x.lumi = beta_mat_male_auto, method = "quantile")

## Perform quantile normalization ...

dim(betaQN_male_auto)

## [1] 446784     27

### male chromosome X
betaQN_male_X <- lumiN(x.lumi = beta_mat_male_X, method = "quantile")

## Perform quantile normalization ...

dim(betaQN_male_X)

## [1] 10618    27

### male chromosome Y
betaQN_male_Y <- lumiN(x.lumi = beta_mat_male_Y, method = "quantile")

## Perform quantile normalization ...

dim(betaQN_male_Y)

## [1] 62 27

4.1.2.2.1 Female plots after normalization

pheno_female_df <- pheno_df %>% filter(sex == "M")
matboxplot(betaQN_female_auto,
           groupFactor = pheno_female_df$sex,
           xaxt = "n",
           main = "Beta Values (autosomes) - after normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_female_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

matboxplot(betaQN_female_X,
           groupFactor = pheno_female_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome X) - after normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_female_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

4.1.2.2.2 Male plots after normalization

pheno_male_df <- pheno_df %>% filter(sex == "M")
matboxplot(betaQN_male_auto,
           groupFactor = pheno_male_df$sex,
           xaxt = "n",
           main = "Beta Values (autosomes) - after normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_male_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

matboxplot(betaQN_male_X,
           groupFactor = pheno_male_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome X) - after normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_male_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

matboxplot(betaQN_male_Y,
           groupFactor = pheno_male_df$sex,
           xaxt = "n",
           main = "Beta Values (chromosome Y) - after normalization")
legend('bottom',
       paste0("sex ", levels(as.factor(pheno_male_df$sex))),
       col = c(1:7), lty = 1, lwd = 3, cex = 0.70)

4.2 BMIQ

library(wateRmelon)
library(RPMM)
library(sesame)
library(sesameData)
library(sm)

4.2.1 test percentage of type I and type II probes then normalize the matrix

4.2.1.1 on autosomes

4.2.1.1.1 for females

### Order annotation in the same order as beta matrix
annotType <- sesameDataGet("HM450.hg19.manifest")
annotType$designTypeNumeric <- ifelse(annotType$designType == "I",1,2)

### Density plot for type I and type II probes
betaQNCompleteCol1_female_auto <- betaQN_female_auto[complete.cases(betaQN_female_auto[,1]), ]
annotTypeCompleteCol1_female_auto <- annotType[row.names(betaQNCompleteCol1_female_auto), ]

sm.density.compare(
    betaQNCompleteCol1_female_auto[,1],
    annotTypeCompleteCol1_female_auto$designTypeNumeric) +
  title(main = "Density plot for type I and type II probes on female autosomes")

## integer(0)

### Summary table
type12 <- annotType$designTypeNumeric[match(rownames(betaQN_female_auto),names(annotType))]
table(type12)

## type12
##      1      2 
## 125484 321300

set.seed(946)
doParallel::registerDoParallel(cores = 4)
betaQN_BMIQ_female_auto <- plyr::aaply(
  betaQN_female_auto, 2,
  function(x){
    norm_ls <- BMIQ(x, design.v = type12, plots = FALSE)
    return (norm_ls$nbeta)
  },.progress = "time",.parallel = TRUE
) %>% t()

## Progress disabled when using parallel plyr

colnames(betaQN_BMIQ_female_auto) <- substr(colnames(betaQN_BMIQ_female_auto),1,stringr::str_length(pheno_df$sample) %>% unique)

4.2.1.1.2 for males

### Density plot for type I and type II probes
betaQNCompleteCol1_male_auto <- betaQN_male_auto[complete.cases(betaQN_male_auto[,1]), ]
annotTypeCompleteCol1_male_auto <- annotType[row.names(betaQNCompleteCol1_male_auto), ]

sm.density.compare(
  betaQNCompleteCol1_male_auto[,1],
  annotTypeCompleteCol1_male_auto$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on male autosomes")

### Summary table
type12 <- annotType$designTypeNumeric[match(rownames(betaQN_male_auto),names(annotType))]
table(type12)

## type12
##      1      2 
## 125484 321300

set.seed(946)
doParallel::registerDoParallel(cores = 4)
betaQN_BMIQ_male_auto <- plyr::aaply(
  betaQN_male_auto, 2,
  function(x){
    norm_ls <- BMIQ(x, design.v = type12, plots = FALSE)
    return (norm_ls$nbeta)
  },.progress = "time",.parallel = TRUE
) %>% t()

## Progress disabled when using parallel plyr

colnames(betaQN_BMIQ_male_auto) <- substr(colnames(betaQN_BMIQ_male_auto),1,stringr::str_length(pheno_df$sample) %>% unique)

4.2.1.2 on chromosome X

4.2.1.2.1 for females

### Density plot for type I and type II probes
betaQNCompleteCol1_female_X <- betaQN_female_X[complete.cases(betaQN_female_X[,1]), ]
annotTypeCompleteCol1_female_X <- annotType[row.names(betaQNCompleteCol1_female_X), ]

sm.density.compare(
  betaQNCompleteCol1_female_X[,1],
  annotTypeCompleteCol1_female_X$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on female chromosome X")

### Summary table
type12 <- annotType$designTypeNumeric[match(rownames(betaQN_female_X),names(annotType))]
table(type12)

## type12
##    1    2 
## 2916 7702

4.2.1.2.2 for males

### Density plot for type I and type II probes
betaQNCompleteCol1_male_X <- betaQN_male_X[complete.cases(betaQN_male_X[,1]), ]
annotTypeCompleteCol1_male_X <- annotType[row.names(betaQNCompleteCol1_male_X), ]

sm.density.compare(
  betaQNCompleteCol1_male_X[,1],
  annotTypeCompleteCol1_male_X$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on male chromosome X")

### Summary table
type12 <- annotType$designTypeNumeric[match(rownames(betaQN_male_X),names(annotType))]
table(type12)

## type12
##    1    2 
## 2916 7702

4.2.1.3 on chromosome Y

### Density plot for type I and type II probes
betaQNCompleteCol1_male_Y <- betaQN_male_Y[complete.cases(betaQN_male_Y[,1]), ]
annotTypeCompleteCol1_male_Y <- annotType[row.names(betaQNCompleteCol1_male_Y), ]

sm.density.compare(
  betaQNCompleteCol1_male_Y[,1],
  annotTypeCompleteCol1_male_Y$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on chromosome Y")

### Summary table
type12 <- annotType$designTypeNumeric[match(rownames(betaQN_male_Y),names(annotType))]
table(type12)

## type12
##  1  2 
## 15 47

4.2.2 save density plots and normalized datasets

### save plots
pdf(paste0(data.dir.probes.normalization, "/densityPlotByProbeType.pdf"))

sm.density.compare(
  betaQNCompleteCol1_female_auto[,1],
  annotTypeCompleteCol1_female_auto$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on female autosomes")

sm.density.compare(
  betaQNCompleteCol1_male_auto[,1],
  annotTypeCompleteCol1_male_auto$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on male autosomes")

sm.density.compare(
  betaQNCompleteCol1_female_X[,1],
  annotTypeCompleteCol1_female_X$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on female chromosome X")

sm.density.compare(
  betaQNCompleteCol1_male_X[,1],
  annotTypeCompleteCol1_male_X$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on male chromosome X")

sm.density.compare(
  betaQNCompleteCol1_male_Y[,1],
  annotTypeCompleteCol1_male_Y$designTypeNumeric
)
title(main = "Density plot for type I and type II probes on chromosome Y")

dev.off()

## quartz_off_screen 
##                 2

### combined normalized matrices
betaQN_BMIQ_female <- rbind(betaQN_BMIQ_female_auto, betaQN_female_X)
betaQN_BMIQ_male <- rbind(betaQN_BMIQ_male_auto, betaQN_male_X, betaQN_male_Y)
save(betaQN_BMIQ_female, betaQN_BMIQ_male, pheno_df, file = paste0(data.dir.probes.normalization, "/GASPARONI_QNBMIQ.rda"))

5 Outliers detection - PCA analysis

• Select most variable probes and perform PCA analysis

Input:

• GASPARONI_QNBMIQ.rds
  
• pheno_df.RDS

Output:

• GASPARONI_PCs_usingBetas.csv,
• PCA plots
• GASPARONI_QNBMIQ_PCfiltered.RDS
• pheno_df.RDS

# plotPCA and OrderDataBySd functions
devtools::source_gist("https://gist.github.com/tiagochst/d3a7b1639acf603916c315d23b1efb3e")

## Sourcing https://gist.githubusercontent.com/tiagochst/d3a7b1639acf603916c315d23b1efb3e/raw/a14424662da343c1301b7b2f03210d28d16ae05c/functions.R

## SHA-1 hash of file is ef6f39dc4e5eddb5ca1c6e5af321e75ff06e9362

### merge in group info
## add center and scale
## compare M values vs. beta values
## with and without center / scale
load(paste0(data.dir.probes.normalization, "/GASPARONI_QNBMIQ.rda"))

5.1 for females

pheno_female_df <- pheno_df %>% filter(sex == "F")
pheno_female_df$stage3 <- ifelse(
  pheno_female_df$stage %in% c(0,1,2), "0-2",
  ifelse(pheno_female_df$stage %in% c(3,4), "3-4", "5-6")
)
betaQN_BMIQ_female <- betaQN_BMIQ_female[ , pheno_female_df$sample]

### transform to m values
mvalue_mat <- log2(betaQN_BMIQ_female / (1 - betaQN_BMIQ_female)) #dim: 457402 30

### order matrix by most variable probes on top
betaOrd_mat <- OrderDataBySd(betaQN_BMIQ_female) #dim: 457402 30
mOrd_mat <- OrderDataBySd(mvalue_mat)  #dim: 457402 30

identical(pheno_female_df$sample, colnames(betaOrd_mat))

## [1] TRUE

identical(pheno_female_df$sample, colnames(mOrd_mat))

## [1] TRUE

pca <- prcomp(t(betaOrd_mat[1:50000, ]),
              center = TRUE,
              scale = TRUE)

d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2])

meanPC1 <- mean (d$PC1)
sdPC1   <- sd (d$PC1)

meanPC2 <- mean (d$PC2)
sdPC2   <- sd (d$PC2)

out3sdPC1_1 <- meanPC1 - 3 * sdPC1
out3sdPC1_2 <- meanPC1 + 3 * sdPC1

out3sdPC2_1 <- meanPC2 - 3 * sdPC2
out3sdPC2_2 <- meanPC2 + 3 * sdPC2

d$outlier_PC1[d$PC1 >= out3sdPC1_1 & d$PC1 <= out3sdPC1_2] <- 0
d$outlier_PC1[d$PC1 < out3sdPC1_1 | d$PC1 > out3sdPC1_2] <- 1

d$outlier_PC2[d$PC2 >= out3sdPC2_1 & d$PC2 <= out3sdPC2_2] <- 0
d$outlier_PC2[d$PC2 < out3sdPC2_1 | d$PC2 > out3sdPC2_2] <- 1

readr::write_csv(d, paste0(data.dir.pca, "/GASPARONI_PCs_usingBetas_female.csv"))

### 2.PCA plot
library(ggplot2)
library(ggrepel)

### Beata values
byStage <- plotPCA(
  dataset = "Gasparoni: beta values",
  expSorted_mat = betaOrd_mat,
  pheno = pheno_female_df,
  group_char = "stage3",
  ntop = 50000,
  center = TRUE,
  scale = TRUE
)

### M values
byStage <- plotPCA(
  dataset = "Gasparoni: M values",
  expSorted_mat = mOrd_mat,
  pheno = pheno_female_df,
  group_char = "stage3",
  ntop = 50000,
  center = TRUE,
  scale = TRUE
)

### Filter samples by PCA, save files
noOutliers <- d[which(d$outlier_PC1 == 0 & d$outlier_PC2 == 0), ]
betaQN_BMIQ_PCfiltered_female <- betaQN_BMIQ_female[, rownames(noOutliers)] #dim: 457402 29
saveRDS(betaQN_BMIQ_PCfiltered_female, paste0(data.dir.pca, "Gasparoni_QNBMIQ_PCfiltered_female.RDS"))

pheno_female_df <- pheno_female_df[pheno_female_df$sample %in% rownames(noOutliers),] #dim: 29 7
saveRDS(pheno_female_df, paste0(data.dir.pca, "/pheno_female_df.RDS"))

5.2 for males

pheno_male_df <- pheno_df %>% filter(sex == "M")
pheno_male_df$stage3 <- ifelse(
  pheno_male_df$stage %in% c(0,1,2), "0-2",
  ifelse(pheno_male_df$stage %in% c(3,4), "3-4", "5-6")
)
betaQN_BMIQ_male <- betaQN_BMIQ_male[ , pheno_male_df$sample]

### transform to m values
mvalue_mat <- log2(betaQN_BMIQ_male / (1 - betaQN_BMIQ_male)) #dim: 457464 27

### order matrix by most variable probes on top
betaOrd_mat <- OrderDataBySd(betaQN_BMIQ_male) #dim: 457464 27
mOrd_mat <- OrderDataBySd(mvalue_mat)  #dim: 457464 27

identical(pheno_male_df$sample, colnames(betaOrd_mat))

## [1] TRUE

identical(pheno_male_df$sample, colnames(mOrd_mat))

## [1] TRUE

pca <- prcomp(t(betaOrd_mat[1:50000, ]),
              center = TRUE,
              scale = TRUE)

d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2])

meanPC1 <- mean (d$PC1)
sdPC1   <- sd (d$PC1)

meanPC2 <- mean (d$PC2)
sdPC2   <- sd (d$PC2)

out3sdPC1_1 <- meanPC1 - 3 * sdPC1
out3sdPC1_2 <- meanPC1 + 3 * sdPC1

out3sdPC2_1 <- meanPC2 - 3 * sdPC2
out3sdPC2_2 <- meanPC2 + 3 * sdPC2

d$outlier_PC1[d$PC1 >= out3sdPC1_1 & d$PC1 <= out3sdPC1_2] <- 0
d$outlier_PC1[d$PC1 < out3sdPC1_1 | d$PC1 > out3sdPC1_2] <- 1

d$outlier_PC2[d$PC2 >= out3sdPC2_1 & d$PC2 <= out3sdPC2_2] <- 0
d$outlier_PC2[d$PC2 < out3sdPC2_1 | d$PC2 > out3sdPC2_2] <- 1

readr::write_csv(d, paste0(data.dir.pca, "/GASPARONI_PCs_usingBetas_male.csv"))

### 2.PCA plot
library(ggplot2)
library(ggrepel)

### Beata values
byStage <- plotPCA(
  dataset = "Gasparoni: beta values",
  expSorted_mat = betaOrd_mat,
  pheno = pheno_male_df,
  group_char = "stage3",
  ntop = 50000,
  center = TRUE,
  scale = TRUE
)

### M values
byStage <- plotPCA(
  dataset = "Gasparoni: M values",
  expSorted_mat = mOrd_mat,
  pheno = pheno_male_df,
  group_char = "stage3",
  ntop = 50000,
  center = TRUE,
  scale = TRUE
)

### Filter samples by PCA, save files
noOutliers <- d[which(d$outlier_PC1 == 0 & d$outlier_PC2 == 0), ]
betaQN_BMIQ_PCfiltered_male <- betaQN_BMIQ_male[, rownames(noOutliers)] #dim: 457464 27
saveRDS(betaQN_BMIQ_PCfiltered_male, paste0(data.dir.pca, "Gasparoni_QNBMIQ_PCfiltered_male.RDS"))

pheno_male_df <- pheno_male_df[pheno_male_df$sample %in% rownames(noOutliers),] #dim: 27 7
saveRDS(pheno_male_df, paste0(data.dir.pca, "/pheno_male_df.RDS"))

6 Summary after QC steps

6.1 Data and metadata

nb.samples.pca <- ncol(betaQN_BMIQ_PCfiltered_female) + ncol(betaQN_BMIQ_PCfiltered_male)
nb.female.samples.pca <- ncol(betaQN_BMIQ_PCfiltered_female)
nb.male.samples.pca <- ncol(betaQN_BMIQ_PCfiltered_male)
dim(betaQN_BMIQ_PCfiltered_female)

## [1] 457402     29

dim(betaQN_BMIQ_PCfiltered_male)

## [1] 457464     27

dim(pheno_female_df)

## [1] 29  7

dim(pheno_male_df)

## [1] 27  7

pheno_female_df %>% 
  DT::datatable(filter = 'top',
                style = "bootstrap",
                extensions = 'Buttons',
                options = list(scrollX = TRUE, 
                               dom = 'Bfrtip',
                               buttons = I('colvis'),
                               keys = TRUE, 
                               pageLength = 10), 
                rownames = FALSE,
                caption = "Female samples metadata")

pheno_male_df %>% 
  DT::datatable(filter = 'top',
                style = "bootstrap",
                extensions = 'Buttons',
                options = list(scrollX = TRUE, 
                               dom = 'Bfrtip',
                               buttons = I('colvis'),
                               keys = TRUE, 
                               pageLength = 10), 
                rownames = FALSE,
                caption = "Male samples metadata")

6.2 Numbers of samples and probes removed in each step

df.samples <- data.frame(
  "Number of samples" =  c(nb.samples, 
                           nb.samples.bc.filtered,
                           nb.samples.with.clinical, 
                           nb.samples.pca),
  "Description" = c("total number of samples",
                    "samples with bisulfate conversion > 88",
                    "samples with clinical data",
                    "Samples after PCA"),
  "Difference" = c("-",
                   nb.samples.bc.filtered - nb.samples ,
                   nb.samples.with.clinical - nb.samples.bc.filtered ,
                   nb.samples.pca - nb.samples.with.clinical)
)    
df.samples           

df.female.samples <- data.frame(
  "Number of samples" =  c(nb.female.samples, 
                           nb.female.samples.bc.filtered,
                           nb.female.samples.with.clinical, 
                           nb.female.samples.pca),
  "Description" = c("Female total number of samples",
                    "Female samples with bisulfate conversion > 88",
                    "Female samples with clinical data",
                    "Female samples after PCA"),
  "Difference" = c("-",
                   nb.female.samples.bc.filtered - nb.female.samples ,
                   nb.female.samples.with.clinical - nb.female.samples.bc.filtered,
                   nb.female.samples.pca - nb.female.samples.with.clinical)
)    
df.female.samples           

df.male.samples <- data.frame(
  "Number of samples" =  c(nb.male.samples, 
                           nb.male.samples.bc.filtered,
                           nb.male.samples.with.clinical, 
                           nb.male.samples.pca),
  "Description" = c("Male total number of samples",
                    "Male samples with bisulfate conversion > 88",
                    "Male samples with clinical data",
                    "Male samples after PCA"),
  "Difference" = c("-",
                   nb.male.samples.bc.filtered - nb.male.samples ,
                   nb.male.samples.with.clinical - nb.male.samples.bc.filtered ,
                   nb.male.samples.pca - nb.male.samples.with.clinical)
)    
df.male.samples           

# Create summary table
df.probes <- data.frame("Number of probes" = c(nb.probes,
                                               nb.probes.detectP, 
                                               nb.probes.detectP.cg,
                                               nb.probes.cg.dmrcate,
                                               nb.probes.dmrcate.chrM),
                        "Description" = c("total number of probes in raw data",
                                          "detection P < 0.01",
                                          "only probes that start with cg",
                                          "DMRcate",
                                          "delete probes on chrM"),
                        "Difference" = c("-",
                                         nb.probes.detectP - nb.probes ,
                                         nb.probes.detectP.cg - nb.probes.detectP,
                                         nb.probes.cg.dmrcate - nb.probes.detectP.cg,
                                         nb.probes.dmrcate.chrM - nb.probes.cg.dmrcate)
)
df.probes

df.chrAuto.probes <- data.frame("Number of probes" = c(nb.chrAuto.probes,
                                                       nb.chrAuto.probes.detectP,
                                                       nb.chrAuto.probes.detectP.cg,
                                                       nb.chrAuto.probes.cg.dmrcate),
                        "Description" = c("autosome total number of probes in raw data",
                                          "autosome detection P < 0.01",
                                          "autosome only probes that start with cg",
                                          "autosome DMRcate"),
                        "Difference" = c("-",
                                         nb.chrAuto.probes.detectP - nb.chrAuto.probes ,
                                         nb.chrAuto.probes.detectP.cg - nb.chrAuto.probes.detectP,
                                         nb.chrAuto.probes.cg.dmrcate - nb.chrAuto.probes.detectP.cg)
)
df.chrAuto.probes

df.chrX.probes <- data.frame("Number of probes" = c(nb.chrX.probes,
                                                    nb.chrX.probes.detectP,
                                                    nb.chrX.probes.detectP.cg,
                                                    nb.chrX.probes.cg.dmrcate),
                        "Description" = c("chromosome X total number of probes in raw data",
                                          "chromosome X detection P < 0.01",
                                          "chromosome X only probes that start with cg",
                                          "chromosome X DMRcate"),
                        "Difference" = c("-",
                                         nb.chrX.probes.detectP - nb.chrX.probes ,
                                         nb.chrX.probes.detectP.cg - nb.chrX.probes.detectP,
                                         nb.chrX.probes.cg.dmrcate - nb.chrX.probes.detectP.cg)
)
df.chrX.probes

df.chrY.probes <- data.frame("Number of probes" = c(nb.chrY.probes,
                                                    nb.chrY.probes.detectP,
                                                    nb.chrY.probes.detectP.cg,
                                                    nb.chrY.probes.cg.dmrcate),
                        "Description" = c("chromosome Y total number of probes in raw data",
                                          "chromosome Y detection P < 0.01",
                                          "chromosome Y only probes that start with cg",
                                          "chromosome Y DMRcate"),
                        "Difference" = c("-",
                                         nb.chrY.probes.detectP - nb.chrY.probes ,
                                         nb.chrY.probes.detectP.cg - nb.chrY.probes.detectP,
                                         nb.chrY.probes.cg.dmrcate - nb.chrY.probes.detectP.cg)
)
df.chrY.probes

7 Session information

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       macOS Catalina 10.15.4      
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/New_York            
##  date     2020-05-06                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package                                       * version  date       lib
##  acepack                                         1.4.1    2016-10-29 [1]
##  affy                                            1.64.0   2019-10-29 [1]
##  affyio                                          1.56.0   2019-10-29 [1]
##  annotate                                        1.64.0   2019-10-29 [1]
##  AnnotationDbi                                 * 1.48.0   2019-10-29 [1]
##  AnnotationFilter                                1.10.0   2019-10-29 [1]
##  AnnotationHub                                 * 2.18.0   2019-10-29 [1]
##  askpass                                         1.1      2019-01-13 [1]
##  assertthat                                      0.2.1    2019-03-21 [1]
##  backports                                       1.1.6    2020-04-05 [1]
##  base64                                          2.0      2016-05-10 [1]
##  base64enc                                       0.1-3    2015-07-28 [1]
##  beanplot                                        1.2      2014-09-19 [1]
##  BiasedUrn                                       1.07     2015-12-28 [1]
##  Biobase                                       * 2.46.0   2019-10-29 [1]
##  BiocFileCache                                 * 1.10.2   2019-11-08 [1]
##  BiocGenerics                                  * 0.32.0   2019-10-29 [1]
##  BiocManager                                     1.30.10  2019-11-16 [1]
##  BiocParallel                                  * 1.20.1   2019-12-21 [1]
##  BiocVersion                                     3.10.1   2019-06-06 [1]
##  biomaRt                                         2.42.1   2020-03-26 [1]
##  Biostrings                                    * 2.54.0   2019-10-29 [1]
##  biovizBase                                      1.34.1   2019-12-04 [1]
##  bit                                             1.1-15.2 2020-02-10 [1]
##  bit64                                           0.9-7    2017-05-08 [1]
##  bitops                                          1.0-6    2013-08-17 [1]
##  blob                                            1.2.1    2020-01-20 [1]
##  BSgenome                                        1.54.0   2019-10-29 [1]
##  bsseq                                           1.22.0   2019-10-29 [1]
##  bumphunter                                    * 1.28.0   2019-10-29 [1]
##  callr                                           3.4.3    2020-03-28 [1]
##  cellranger                                      1.1.0    2016-07-27 [1]
##  checkmate                                       2.0.0    2020-02-06 [1]
##  cli                                             2.0.2    2020-02-28 [1]
##  cluster                                       * 2.1.0    2019-06-19 [1]
##  codetools                                       0.2-16   2018-12-24 [1]
##  colorspace                                      1.4-1    2019-03-18 [1]
##  crayon                                          1.3.4    2017-09-16 [1]
##  crosstalk                                       1.1.0.1  2020-03-13 [1]
##  curl                                            4.3      2019-12-02 [1]
##  data.table                                      1.12.8   2019-12-09 [1]
##  DBI                                             1.1.0    2019-12-15 [1]
##  dbplyr                                        * 1.4.3    2020-04-19 [1]
##  DelayedArray                                  * 0.12.3   2020-04-09 [1]
##  DelayedMatrixStats                              1.8.0    2019-10-29 [1]
##  desc                                            1.2.0    2018-05-01 [1]
##  devtools                                        2.3.0    2020-04-10 [1]
##  dichromat                                       2.0-0    2013-01-24 [1]
##  digest                                          0.6.25   2020-02-23 [1]
##  DMRcate                                       * 2.0.7    2020-01-10 [1]
##  DMRcatedata                                   * 2.2.1    2020-02-27 [1]
##  DNAcopy                                         1.60.0   2019-10-29 [1]
##  doParallel                                      1.0.15   2019-08-02 [1]
##  doRNG                                           1.8.2    2020-01-27 [1]
##  dplyr                                         * 0.8.5    2020-03-07 [1]
##  DSS                                             2.34.0   2019-10-29 [1]
##  DT                                              0.13     2020-03-23 [1]
##  edgeR                                           3.28.1   2020-02-26 [1]
##  ellipsis                                        0.3.0    2019-09-20 [1]
##  ensembldb                                       2.10.2   2019-11-20 [1]
##  evaluate                                        0.14     2019-05-28 [1]
##  ExperimentHub                                 * 1.12.0   2019-10-29 [1]
##  fansi                                           0.4.1    2020-01-08 [1]
##  farver                                          2.0.3    2020-01-16 [1]
##  fastmap                                         1.0.1    2019-10-08 [1]
##  FDb.InfiniumMethylation.hg19                  * 2.2.0    2020-03-18 [1]
##  foreach                                       * 1.5.0    2020-03-30 [1]
##  foreign                                         0.8-76   2020-03-03 [1]
##  Formula                                         1.2-3    2018-05-03 [1]
##  fs                                              1.4.1    2020-04-04 [1]
##  genefilter                                      1.68.0   2019-10-29 [1]
##  GenomeInfoDb                                  * 1.22.1   2020-03-27 [1]
##  GenomeInfoDbData                                1.2.2    2020-03-18 [1]
##  GenomicAlignments                               1.22.1   2019-11-12 [1]
##  GenomicFeatures                               * 1.38.2   2020-02-15 [1]
##  GenomicRanges                                 * 1.38.0   2019-10-29 [1]
##  GEOquery                                        2.54.1   2019-11-18 [1]
##  ggplot2                                       * 3.3.0    2020-03-05 [1]
##  ggpubr                                          0.2.5    2020-02-13 [1]
##  ggrepel                                       * 0.8.2    2020-03-08 [1]
##  ggsignif                                        0.6.0    2019-08-08 [1]
##  glue                                            1.4.0    2020-04-03 [1]
##  GO.db                                           3.10.0   2020-03-18 [1]
##  gridExtra                                       2.3      2017-09-09 [1]
##  gtable                                          0.3.0    2019-03-25 [1]
##  gtools                                          3.8.2    2020-03-31 [1]
##  Gviz                                            1.30.3   2020-02-17 [1]
##  HDF5Array                                       1.14.4   2020-04-13 [1]
##  Hmisc                                           4.4-0    2020-03-23 [1]
##  hms                                             0.5.3    2020-01-08 [1]
##  htmlTable                                       1.13.3   2019-12-04 [1]
##  htmltools                                       0.4.0    2019-10-04 [1]
##  htmlwidgets                                     1.5.1    2019-10-08 [1]
##  httpuv                                          1.5.2    2019-09-11 [1]
##  httr                                            1.4.1    2019-08-05 [1]
##  IlluminaHumanMethylation450kanno.ilmn12.hg19  * 0.6.0    2020-03-18 [1]
##  IlluminaHumanMethylation450kmanifest          * 0.4.0    2020-03-18 [1]
##  IlluminaHumanMethylationEPICanno.ilm10b4.hg19   0.6.0    2020-03-18 [1]
##  IlluminaHumanMethylationEPICmanifest            0.3.0    2020-03-18 [1]
##  illuminaio                                    * 0.28.0   2019-10-29 [1]
##  interactiveDisplayBase                          1.24.0   2019-10-29 [1]
##  IRanges                                       * 2.20.2   2020-01-13 [1]
##  iterators                                     * 1.0.12   2019-07-26 [1]
##  jpeg                                            0.1-8.1  2019-10-24 [1]
##  jsonlite                                        1.6.1    2020-02-02 [1]
##  KernSmooth                                      2.23-17  2020-04-26 [1]
##  knitr                                           1.28     2020-02-06 [1]
##  labeling                                        0.3      2014-08-23 [1]
##  later                                           1.0.0    2019-10-04 [1]
##  lattice                                         0.20-41  2020-04-02 [1]
##  latticeExtra                                    0.6-29   2019-12-19 [1]
##  lazyeval                                        0.2.2    2019-03-15 [1]
##  lifecycle                                       0.2.0    2020-03-06 [1]
##  limma                                         * 3.42.2   2020-02-03 [1]
##  locfit                                        * 1.5-9.4  2020-03-25 [1]
##  lumi                                          * 2.38.0   2019-10-29 [1]
##  magrittr                                        1.5      2014-11-22 [1]
##  MASS                                            7.3-51.6 2020-04-26 [1]
##  Matrix                                          1.2-18   2019-11-27 [1]
##  matrixStats                                   * 0.56.0   2020-03-13 [1]
##  mclust                                          5.4.6    2020-04-11 [1]
##  memoise                                         1.1.0    2017-04-21 [1]
##  methylumi                                     * 2.32.0   2019-10-29 [1]
##  mgcv                                            1.8-31   2019-11-09 [1]
##  mime                                            0.9      2020-02-04 [1]
##  minfi                                         * 1.32.0   2019-10-29 [1]
##  missMethyl                                      1.20.4   2020-01-28 [1]
##  multtest                                        2.42.0   2019-10-29 [1]
##  munsell                                         0.5.0    2018-06-12 [1]
##  nleqslv                                         3.3.2    2018-05-17 [1]
##  nlme                                            3.1-147  2020-04-13 [1]
##  nnet                                            7.3-14   2020-04-26 [1]
##  nor1mix                                         1.3-0    2019-06-13 [1]
##  openssl                                         1.4.1    2019-07-18 [1]
##  org.Hs.eg.db                                  * 3.10.0   2020-03-18 [1]
##  permute                                         0.9-5    2019-03-12 [1]
##  pillar                                          1.4.3    2019-12-20 [1]
##  pkgbuild                                        1.0.7    2020-04-25 [1]
##  pkgconfig                                       2.0.3    2019-09-22 [1]
##  pkgload                                         1.0.2    2018-10-29 [1]
##  plyr                                            1.8.6    2020-03-03 [1]
##  png                                             0.1-7    2013-12-03 [1]
##  preprocessCore                                  1.48.0   2019-10-29 [1]
##  prettyunits                                     1.1.1    2020-01-24 [1]
##  processx                                        3.4.2    2020-02-09 [1]
##  progress                                        1.2.2    2019-05-16 [1]
##  promises                                        1.1.0    2019-10-04 [1]
##  ProtGenerics                                    1.18.0   2019-10-29 [1]
##  ps                                              1.3.2    2020-02-13 [1]
##  purrr                                           0.3.4    2020-04-17 [1]
##  quadprog                                        1.5-8    2019-11-20 [1]
##  quantro                                       * 1.20.0   2019-10-29 [1]
##  R.methodsS3                                     1.8.0    2020-02-14 [1]
##  R.oo                                            1.23.0   2019-11-03 [1]
##  R.utils                                         2.9.2    2019-12-08 [1]
##  R6                                              2.4.1    2019-11-12 [1]
##  randomForest                                    4.6-14   2018-03-25 [1]
##  rappdirs                                        0.3.1    2016-03-28 [1]
##  RColorBrewer                                    1.1-2    2014-12-07 [1]
##  Rcpp                                            1.0.4.6  2020-04-09 [1]
##  RCurl                                           1.98-1.2 2020-04-18 [1]
##  readr                                           1.3.1    2018-12-21 [1]
##  readxl                                          1.3.1    2019-03-13 [1]
##  remotes                                         2.1.1    2020-02-15 [1]
##  reshape                                         0.8.8    2018-10-23 [1]
##  reshape2                                      * 1.4.4    2020-04-09 [1]
##  rhdf5                                           2.30.1   2019-11-26 [1]
##  Rhdf5lib                                        1.8.0    2019-10-29 [1]
##  rlang                                           0.4.5    2020-03-01 [1]
##  rmarkdown                                       2.1      2020-01-20 [1]
##  rngtools                                        1.5      2020-01-23 [1]
##  ROC                                           * 1.62.0   2019-10-29 [1]
##  rpart                                           4.1-15   2019-04-12 [1]
##  RPMM                                          * 1.25     2017-02-28 [1]
##  rprojroot                                       1.3-2    2018-01-03 [1]
##  Rsamtools                                       2.2.3    2020-02-23 [1]
##  RSQLite                                         2.2.0    2020-01-07 [1]
##  rstudioapi                                      0.11     2020-02-07 [1]
##  rtracklayer                                     1.46.0   2019-10-29 [1]
##  ruv                                             0.9.7.1  2019-08-30 [1]
##  S4Vectors                                     * 0.24.4   2020-04-09 [1]
##  scales                                        * 1.1.0    2019-11-18 [1]
##  scrime                                          1.3.5    2018-12-01 [1]
##  sesame                                        * 1.4.0    2019-10-29 [1]
##  sesameData                                    * 1.4.0    2019-11-05 [1]
##  sessioninfo                                     1.1.1    2018-11-05 [1]
##  shiny                                           1.4.0.2  2020-03-13 [1]
##  siggenes                                        1.60.0   2019-10-29 [1]
##  sm                                            * 2.2-5.6  2018-09-27 [1]
##  statmod                                         1.4.34   2020-02-17 [1]
##  stringi                                         1.4.6    2020-02-17 [1]
##  stringr                                         1.4.0    2019-02-10 [1]
##  SummarizedExperiment                          * 1.16.1   2019-12-19 [1]
##  survival                                        3.1-12   2020-04-10 [1]
##  testthat                                        2.3.2    2020-03-02 [1]
##  tibble                                          3.0.1    2020-04-20 [1]
##  tidyr                                           1.0.2    2020-01-24 [1]
##  tidyselect                                      1.0.0    2020-01-27 [1]
##  TxDb.Hsapiens.UCSC.hg19.knownGene             * 3.2.2    2020-03-18 [1]
##  usethis                                         1.6.0    2020-04-09 [1]
##  VariantAnnotation                               1.32.0   2019-10-29 [1]
##  vctrs                                           0.2.4    2020-03-10 [1]
##  wateRmelon                                    * 1.30.0   2019-10-29 [1]
##  wheatmap                                        0.1.0    2018-03-15 [1]
##  withr                                           2.2.0    2020-04-20 [1]
##  xfun                                            0.13     2020-04-13 [1]
##  XML                                             3.99-0.3 2020-01-20 [1]
##  xml2                                            1.3.2    2020-04-23 [1]
##  xtable                                          1.8-4    2019-04-21 [1]
##  XVector                                       * 0.26.0   2019-10-29 [1]
##  yaml                                            2.2.1    2020-02-01 [1]
##  zlibbioc                                        1.32.0   2019-10-29 [1]
##  source        
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.2)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  CRAN (R 3.6.2)
##  CRAN (R 3.6.0)
##  Bioconductor  
##  CRAN (R 3.6.0)
##  Bioconductor  
## 
## [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library